The joint model for the hazard functions for recurrent event $r_i(.)$ and death $\lambda_i(.)$ is :
$$r_i(t|\omega_i)=\omega_ir_0(t)exp(\bold{\beta_1^{'}Z_i(t)})$$ $$\lambda_i(t|\omega_i)=\omega_i^{\alpha}\lambda_0(t)exp(\bold{\beta_2^{'}Z_i(t)})$$
where $r_0(t)$ (resp. $\lambda_0(t)$) is the recurrent (resp. terminal) event baseline hazard function, $\bold{\beta_1}$ (resp. $\bold{\beta_2}$) the regression coefficient vector, $\bold{Z_i(t)}$ the covariate vector. And where $\omega_i\sim\bold{\Gamma}(\frac{1}{\theta},\frac{1}{\theta})$ is iid.
INITIAL VALUES
When a joint frailty model is fitted, at first, the splines coefficients and the regression coefficients are initialized at 0.5, and the variance of the frailty term $\theta$ and the coefficient $\alpha$ associated in the death hazard function are initialized at 1. Thus, the next step is the maximisation of the penalized log-likelihood.
PARAMETERS LIMIT VALUES
As frailtypack is written in Fortran 77 some parameters had to be hard coded in.The default values of these parameters are, with the corresponding variable namein the fortran code between brackets.
maximum number of observations(ndatemax): 30000 maximum number of groups(ngmax): 1000 maximum number of subjects(nsujetmax): 15000 maximum number of parameters (npmax) :50 maximum number of covariates (nvarmax):50 If these parameters are not large enough (an error message will let you know this), you need to reset them in joint.f and recompile.
V. Rondeau, D Commenges, and P. Joly (2003). Maximum penalized likelihood estimation in a gamma-frailty model. Lifetime Data Analysis 9, 139-153.
D. Marquardt (1963). An algorithm for least-squares estimation of nonlinear parameters. SIAM Journal of Applied Mathematics, 431-441.
summary.jointPenal,
print.jointPenal,
plot.jointPenal,
dataJoint,
terminal,
cluster### Joint model (recurrent and terminal events) with 2 covariates ###
### on a simulated dataset ###
data(dataJoint)
modJoint<-frailtyPenal(Surv(time.entry,time.end,status)~cluster(id)+var1+var2
+terminal(status.terminal),
formula.terminalEvent=~var1,
data=dataJoint,n.knots=7,
kappa1=1, kappa2=1, Frailty = TRUE, joint=TRUE,
recurrentAG=TRUE)
print(modJoint)
summary(modJoint)
plot(modJoint)
# It takes around 1 minute to converge #Run the code above in your browser using DataLab